Interpreting the results of searches for gravitational waves from coalescing binaries 
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We introduce a method based on the loudest event statistic to calculate an upper limit or interval on the 
astrophysical rate of binary coalescence. The calculation depends upon the sensitivity and noise background 
of the detectors, and a model for the astrophysical distribution of coalescing binaries. There are significant 
uncertainties in the calculation of the rate due to both astrophysical and instrumental uncertainties as well as 
errors introduced by using the post-Newtonian waveform to approximate the full signal. We catalog these 
uncertainties in detail and describe a method for marginalizing over them. Throughout, we provide an example 
based on the initial LIGO detectors. 



I. INTRODUCTION 

The first generation of gravitational wave interferometric 
detectors have achieved, or are approaching, their design sen- 
sitivities mS.!!!]. One of the most promising sources of grav- 
itational waves for these detectors are those emitted during 
the coalescence of a binary system composed of neutron stars 
or black holes. The initial detectors will be sensitive to binary 
neutron star (BNS) coalescences as far away as the Virgo clus- 
ter, while for binary black holes (BBH) the reach is as great 
as 100 Mpc. Thus, with a year of data, there is a possibil- 
ity of detecting gravitational waves from these sources. Even 
in the absence of a detection, the upper limits which will be 
placed on the rates of binary coalescences will begin to be 
astrophysically interesting. Following the initial searches, the 
detectors will be upgraded to enhanced and advanced configu- 
rations where the sensitivity will increase by factors of around 
two and ten respectively over the initial detectors. In this pa- 
per we describe a method which can be used for providing an 
astrophysical interpretation of the results of a search for com- 
pact binary coalescence. In the absence of a detection, this 
will result in an upper limit on the rate of such events. If a 
signal is observed, then a rate interval will be calculated. The 
rate limit calculation makes use of the loudest event statistic, 
first described in |41j and elaborated upon in ISi] . 

The rate of gravitational waves observable at the detectors 
depends upon the detector sensitivity and also the astrophysi- 
cal model for the source in question. The standard assumption 
10] is that the rate of binary coalescence in a galaxy follows 
the blue light luminosity, which depends upon the star for- 
mation rate. At cosmological distances, it is appropriate to 
treat the universe as homogeneous and isotropic. However, 
for the initial and enhanced gravitational-wave detectors, lo- 
cal anisotropics can have a significant effect. In particular, for 
initial LIGO and Virgo, a significant fraction of the available 
luminosity comes from the Virgo cluster ff]. Therefore, it is 
critically important to generate an accurate distribution of the 
blue light in the local universe. In Ref. J^l a catalog of local 
galaxies has been constructed for precisely this purpose. In 
this paper, we describe in detail how this galaxy catalog can be 
used to calculate the total luminosity to which a gravitational- 
wave detector is sensitive. In addition to the galaxy distribu- 



tion, this is dependent upon evaluating the ability of a search 
to detect a binary with given parameters. We describe how 
this search efficiency can be calculated and folded together 
with the galaxy distribution to obtain a total luminosity for a 
given search. 

The amplitude of the gravitational-wave signal emitted by 
a coalescing binary is dependent upon the component masses 
of the binary. Thus, a rate calculation (upper limit or inter- 
val) obtained from a search of gravitational-wave data will 
be sensitive to the astrophysical model for the distribution of 
the component masses of the binary. Although several neu- 
tron stars in binaries, and many isolated neutron stars, have 
been observed as pulsars, there is still significant uncertainty 
in the mass distribution i^. Furthermore, given the lack of 
observed neutron star-black hole or binary black hole sys- 
tems, the distribution of component masses of these systems 
is highly uncertain. With this in mind, we describe a simple 
method whereby mass dependent rates can be obtained. 

In order to claim the detection of gravitational waves, it is 
vitally important to have a good understanding of the distribu- 
tion of events which are due to instrumental noise. The output 
from a search for gravitational-wave transients is a list of can- 
didate events which survive all thresholds applied during the 
search. These candidate events can contain both background, 
noise events as well as true gravitational- wave signals. With- 
out a good understanding of the noise background, it will be 
impossible to determine that an event is due to a gravitational- 
wave signal. We shall see in detail later that a good estimate 
of the noise background is equally important for calculating 
rate limits or intervals. 

There are numerous uncertainties affecting our understand- 
ing of the astrophysics, instrumental sensitivities and back- 
ground estimates required in interpreting the results^ The un- 
certainties in the galaxy distribution are detailed in |8i]. While 
the sky position of nearby galaxies is well known, their dis- 
tance from the earth and blue light luminosity is measured 
less accurately. The search for gravitational waves from coa- 
lescing binaries makes use of waveform templates calculated 
using the post-Newtonian approximation to general relativ- 
ity. This will lead to some discrepancy between the waveform 
used in the search and true physical waveform, particularly 
close to coalescence. In addition, the spin of the binaries is 
neglected in current searches ifioll . However, neutron stars 



2 



and black holes occurring in binary systems are expected to 
be spinning. The background of noise events is estimated by 
time shifting the data from each instrument relative to the oth- 
ers. This provides a reasonable, but not perfect estimate of the 
background rate. Furthermore, instrumental calibration uncer- 
tainties will affect the reported sensitivity of the instrument. 
All of these will have an effect on the calculated rates. We 
provide a detailed discussion of these systematic uncertainties 
and describe a method by which these can be incorporated 
into the final rate statements. 

The layout of the paper is as follows. First, in Section Ull 
we briefly review the techniques used in searching for grav- 
itational waves from coalescing binaries. In Section Hill we 
discuss the loudest event statistic, and demonstrate its applica- 
tion to a search for binary inspiral. In Section|IV]we describe 
the various systematic effects which can effect the rate upper 
limit, and describe a method for marginalizing over these un- 
certainties. Finally, in SectionlVlwe summarize the results. 

Throughout the paper, the methods discussed are applied 
to an illustrative example. This example involves simulated 
results from two detectors operating at the initial LIGO sen- 
sitivity for one year More concretely we choose instruments 
with a binary neutron star horizon distance (the distance at 
which an optimally oriented and located system would pro- 
duce a signal to noise ratio (SNR) of 8 in a single detector) of 
35 Mpc. This is consistent with sensitivities achieved by the 
Hanford and Livingston 4km detectors during the latter parts 
of the S5 science run. Furthermore, we assume that the detec- 
tors and search are such that the SNR threshold can be set at 
5.5 and noise background will produce an expected rate of one 
event per year at a combined SNR of 10, or equivalently about 
7 in each instrument. This value is somewhat greater than the 
combined SNR of 8 often assumed (see for example ifTTll ) and 
reflects a realistic value given the non-stationarity of the initial 
detectors. Although the example presented here involves sim- 
ulated results, the methods are easily generalized to searches 
of real data. Indeed, the reported upper limits from search- 
ing the data from the third and fourth LIGO science runs ifioll 
were obtained using the methods described here. 



II. DETAILS OF INSPIRAL SEARCHES 



The gravitational waveform emitted by a coalescing binary 
can be calculated to high accuracy within the post-Newtonian 
expansion of general relativity. In Section HTaI we sketch the 
details of the inspiral waveform, paying particular attention to 
those quantities which have a direct impact on the detectabil- 
ity of the waveform at a given detector. Then, in Section HlBl 
we describe the matched-filter search method which is used in 
searching for the waveform. We finish with a brief description 
of the extension to multi-detector searches. 



A. The waveform 

The gravitational- wave strain incident on an interferometer is 
given by 

h{t) ^ F+{0,cf),^b)h+{t) + F^{0,^,^)h^{t) , (1) 

where and Fx are the detector response functions and h+, 
hx are the two polarizations of the gravitational wave. The de- 
tector response functions depend upon the sky location {9, cj)) 
and polarization angle tp of the source relative to the detector 
according to llT3l 

F+ = — -(1 + cos^ 61) cos2(/)cos2?/' — cos6'sin2(^sin2'0 , 
Fx = -(1 + cos^ 6') cos 20 sin 2'(/; — cos 6* sin 20 cos 2^/1.(2) 

For a low-mass binary coalescence, the portion of the wave- 
form which is available to the LIGO and Virgo detectors can 
be calculated using the restricted post-Newtonian expansion 
to a high accuracy |ll [M [II El ll2l [III, JJ. 20, 22, 21 . 
For restricted waveforms, only the leading order term in the 
amplitude is used, while post-Newtonian corrections to the 
phase are retained. The waveform depends upon the masses 
and spins of the binary components, the orientation and dis- 
tance of the binary relative to the detector. In this paper, we 
restrict attention to searches which make use of waveforms 
appropriate for binaries without spin. In Section lTV BI we esti- 
mate the effect of using non-spinning templates in the search 
for spinning waveforms. 

The post-Newtonian expression for the binary inspiral 
waveform can be substituted into Eq. ([T]i above to obtain an 
alternative expression for the gravitational-wave strain at the 
detector ll24il25ll . In this form, the restricted waveform is ex- 
pressed in terms of the two phases — he, the cosine phase, 
and hs, the sine phase — of the chirp waveform, and an over- 
all amplitude and phase factor. By doing so, there is a clear 
split between the overall scaling factors which depend upon 
the distance, sky location and orientation of the binary and the 
mass dependent time evolution of the waveform. The wave- 
form is written 

h{t) ^^^^[hc{t) cos ^ + h,{t) sin <^>] , (3) 

where the amplitude factor Dcff is known as the effective dis- 
tance to the binary and <1> is the coalescence phase as observed 
at the detector. Both the effective distance and coalescence 
phase are determined by the location and orientation of the 
binary system relative to the detector More specifically, the 
effective distance is defined as ll26ll 

DoS = , , (4) 

JFlil + cos2 t)2/4 + F2 cos2 i 

where r is the physical distance to the binary, t is the inclina- 
tion angle of the binary system and F+, Fx are given in (|2|i. 
Similarly, the phase angle $ is dependent upon the sky loca- 
tion, polarization, inclination and also the coalescence phase 
(po of the binary. 
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Define the Fourier transform of h{t) by 

/•OO 

h{f)= / h{t)e-^^'f'dt . 



(5) 



The sine and cosine phases of the binary inspiral waveform 
are dependent upon the component masses. In the frequency 
domain, they are obtained using the stationary phase approxi- 
mation to the post-Newtonian expansion ll27ll . In this approx- 
imation, hs{f) = ihc{f), and 



GM 



5/6 



e(/ - /max)/-'/'e^*(^^-^'"' , (6) 



where M is an overall (known) scaling, and we have intro- 
duced the chirp mass M. = Alrj^^^, where M = {nii + 7712) 
is the total mass and 77 = mim2/M^ is the symmetric mass 
ratio. The phase evolution is governed by ^'(/; A4, rf) which 
also depends upon the masses of the binary system. The wave- 
form terminates at a frequency /max- Physically, the post- 
Newtonian waveform is expected to break down around the 
frequency where the evolution changes from a slow inspiral 
to a rapid merger A reasonable estimate of this frequency 
is given by the innermost stable circular orbit (ISCO) of the 
Schwarzschild spacetime with the same total mass M, 



/isco = ^ 1600Hz. 



(7) 



In principle, the waveform observed at the detector for a 
non-spinning binary system depends upon eight parameters: 
the masses of the two binary components, the physical dis- 
tance r, the sky location and polarization {9, (/>, 1/'), the incli- 
nation angle t and the coalescence phase 00 . However, it is 
clear from Eq. (|6]l that the six quantities describing the loca- 
tion and orientation of the binary affect the strain observed 
at a single detector only in the combinations Doff and $. Fur- 
thermore, the ability to distinguish a gravitational wave from a 
coalescing binary above the background noise is independent 
of the coalescence phase at the detector Finally, it is only the 
chirp mass which affects the amplitude of the waveform (the 
mass ratio rj does affect the phase evolution). Therefore, the 
detection efficiency will depend primarily upon the effective 
distance Dcff and the chirp mass A4. This observation will be 
used to greatly simplify the rate calculation. 



B. Inspiral Search Method 

Since the inspiral waveform is well known, the standard 
matched filtering technique is used to distinguish signal from 
noise in a single detector ll28ll . Here, we provide a very brief 
review of the search implementation, further details are avail- 
able in Ref. [24, 25]. The gravitational waveform from a co- 
alescing binary given in Eq. (|6]l depends upon the masses, 
effective distance and coalescence phase of the binary. The 
two mass dimensions are searched by covering the mass space 
with a template bank which guarantees that for any signal 
there is a good overlap between the waveform and the closest 



template 112911 . As discussed below, the coalescence phase are 
analytically maximized over in the matched filtering process 
and the distance is measured. 
The output of a detector is 



5(t) = n{t) + h{t) 



(8) 



where n{t) is the instrumental noise and h{t) is some, possi- 
bly absent, signal. The sensitivity of the instrument is charac- 
terized by the (one-sided) power spectrum S{f), given as 



(9) 



where the tilde represents the Fourier transform, and the angle 
brackets denote the expectation value over the noise. In order 
to construct the matched filter, we introduce an inner product 



{a\b) := 4Rc / df 



S{\f\) 



(10) 



The sensitivity of the detector to a given signal is character- 
ized by 



ihc\hc) , 



(11) 



where depends upon the noise curve of the instrument as 
well as the masses of the binary components (recall that he is 
the waveform at an effective distance of 1 Mpc). The SNR is 
defined as 



2 

P = 



{S\h,f + {S\hsf 



(12) 



The analytic maximization over the unknown phase angle has 
already been performed by including the outputs from the two 
orthogonal filters, he and hs, while the amplitude of the signal 
determines the value of the SNR. 

In the absence of a signal, the SNR squared is distributed 
with two degrees of freedom — one for each of the filters. 
Thus for a single trial, the probability of obtaining an SNR 
greater than is 



2^ _ „-p:/2 



(13) 



If the detector output contains a signal h{t), we characterize 
its amplitude by 



pi := {h\h) = 



(14) 



In this case, the SNR squared is distributed as a non-central 
distribution with two degrees of freedom and a non-centrality 
parameter p\. Thus, the expected SNR squared is + 2 while 
the variance is 4(p^ + 1). For SNRs well above unity, the 
expected SNR is close to ph. 

In the course of a gravitational wave search, we calculate 
the SNR for every mass template in the bank and for every 
time point. This is used to construct a list of times and asso- 
ciated mass parameters at which the SNR exceeds some pre- 
determined threshold. These candidate events may be due to 
instrumental noise or gravitational waves. In searching a year 
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of data over a wide range of masses, we obtain a background 
of events due to noise with a distribution consistent with 
(fTsl l. In addition, the data from the detectors contains non- 
stationarities which also produce events with high SNR. The 
ability to reduce this background of non-gravitational wave 
induced events affects the sensitivity of the search. There are 
many techniques employed to achieve this jsoll. The most 
powerful tool is a consistency test between detectors — the 
arrival time of the signal should be consistent, within the 
light travel time between the sites, and the mass parameters 
should agree within the search accuracy. In addition, sig- 
nal consistency tests, such as the 113 111 and tests 113211 
are utilized, and an "effective SNR", constructed using the 
SNR and information, is used to distinguish signal from 
noise ifioilsoll . After applying these tests, the typical loudest 
surviving background events for BNS searches have a com- 
bined SNR p = \/ p\ of around 10. For higher mass 
systems, the waveforms spend a shorter time in the sensitive 
band. Consequently, it is more difficult to distinguish them 
from noise non-stationarities whence the background extends 
to higher SNR. 



than p is' 

P^(p)=.e-«^-('')^. (15) 

where Cl{p) is the total luminosity to which the search is 
sensitive and T is the duration of the search. Similarly, the 
probability of obtaining zero accidental noise (or background) 
events with an SNR greater than p is denoted Pb{p)- There- 
fore, the overall probability of obtaining no events with SNR 
greater than p is 

P{p\B,R,T)^PB{p)e-~'^'^^^p'^'^ . (16) 

Given that no event was observed with SNR greater than 
that of the loudest event pm, a straightforward application of 
Bayes theorem leads to the posterior probability distribution 
for the rate T, B) which depends upon the loudest 

event, amount of time searched and the accidental rate (en- 
coded as 'B') jsl]. In addition, it will depend upon the prior 
probability distribution for the rate, denoted p{R). In the 
absence of previous knowledge of the rate, a uniform prior 
p{R) = const, is appropriate, while if a previous search has 
been performed, the posterior of that search is naturally used 
as the prior for the next search. The expression for the poste- 
rior distribution is 



m. RATE CALCULATIONS FOR INSPIRAL SEARCHES 



Let us assume that a search for coalescing binaries has been 
performed on a stretch of data from gravitational-wave detec- 
tors. We would like to use the search results to make a state- 
ment about the rate of binary coalescences in the universe. 
This can be done by making use of the loudest event statistic, 
as described in ||3,|5|]. The result will depend upon the astro- 
physical model for the distribution of binary coalescences. To 
roceed, we make the simple, yet astrophysically reasonable 
assumption that binary coalescences occur only in galax- 
ies and furthermore the rate of binary coalescence in a given 
galaxy is directly proportional to the blue light luminosity of 
that galaxy. This assumption is justified by the fact that both 
the star formation rate and supernova rate are proportional to 
the blue light luminosity. The result from a search for coalesc- 
ing binaries, in the absence of a detection, will be a bound on 
the rate R of binary inspirals per year per Lio = 10'°Lb 0, 
where Lb.q is the blue light luminosity of the sun. Recent 
papers have suggested that due to the large delay between for- 
mation and coalescence for binaries, this simple assumption 
will underestimate the contribution from elliptical galaxies, 
particularly for BBH rates Issll . The formalism presented be- 
low can be modified in a straightforward manner to take this 
into account. 

The loudest event statistic makes use of the fact that no 
events with an SNR greater than that of the loudest event p„-^ 
occurred in the search. For an inspiral rate R, the probability 
of there being no gravitational-wave signals with SNR greater 



p(i?|p™,T,B) ocp(i?) 



l+KRCL{Pm)T 



1 + A 



(17) 



Here, A is a measure of the likelihood that the loudest event 
is a due to a gravitational wave, rather than from instrumental 
noise. The expression for A is 



A = 



\C'L{p.n)\ 
P'BiPrn) 



ChiPm) 



PB{Pin) 



(18) 



where the derivatives are taken with respect to p. The likeli- 
hood depends upon the population through the total luminos- 
ity CL{pm) and its derivative with respect to SNR, C£(p„i). 
Similarly it depends upon the background distribution of noise 
events through PB{pm) and its derivative. For practical cal- 
culational purposes, it is often useful to write 



riB 



where np ~ 



\c'M\ 

ChiPm) 



Ub 



Pb{p^) 

PsiPm) 



(19) 

By doing this, the dependence upon the estimated background 
is confined to ub while the quantity np depends only upon the 
astrophysical population model. 

The posterior distribution obtained in Eq. ( fTTI ) can be used 
to calculate a Bayesian upper limit on the rate. The upper limit 
i?* for a given confidence level (a) is obtained by evaluating 



p{R\prn,T, B) . 



(20) 



' The notation F is used to denote "foreground" or gravitational-wave 
events, in contrast to background B events associated to instrumental noise. 
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Assuming a uniform prior on the rate, the upper limit is given 



by 



1 — a = e' 



A 



1 + A 



R*TCl{p 



(21) 

In the case where the loudest event candidate is most likely 
due to the background A ^ and the upper limit becomes 
i?90% = 2.3/[rCL(p™)]. 

In the limit that A ^ oo, the event is almost definitely due 
to a gravitational wave rather than from the noise background. 
In this limit, the probability distribution in Eq. ( [TtI i is peaked 
away from zero, whence it makes sense to construct a rate 
interval (as described in ^) rather than a rate upper limit. 
This is done by performing the integral in ( l20b from i?inin to 
i?,nax- Alternatively, one could choose to still place an upper 
hmit, in which case i?9o% = 3.9/[TCL(pm)]. 

In the remainder of this section, we discuss how the quanti- 
ties A and C^, can be calculated. In Section lTlI Al we describe 
the estimation of the noise background and the evaluation of 
ub- Then, in Section lTlI Bl we describe the calculation of the 
luminosity Cl and consequently iip- Finally, in Section lTlI CI 
we combine these results to obtain an expression for the upper 
limit. 



A. Estimating the Noise Background 

In order to calculate an upper limit, we require an estimate 
of the background of events due to noise in the detectors. For 
perfectly Gaussian, stationary detectors, this can be calculated 
theoretically using the known distribution for the SNR. How- 
ever, real detectors suffer from non-stationarities and time 
varying sensitivity. Thus, for a search involving several detec- 
tors, the background is most readily estimated by time shifting 
the data from the detectors relative to one another, and then 
searching for events which are coincident in time and mass 
between the detectors. If the time shifts are greater than the 
light travel time and length of the signal, then the time shifted 
coincidences cannot be due to gravitational waves, and are 
therefore expected to give a good estimate of the background. 
Each time shift will have a loudest event, which we assume 
to be drawn from the actual background distribution for the 
loudest event, pb{p) P'siP)- Therefore, by performing 
a large number of time shifts, we obtain a good sampling of 
Pb{p) from which it is straightforward to obtain the cumula- 
tive distribution Pb{p) and ns = PsiPm) / PsiPm)- 

As an example, consider a pair of detectors whose noise 
output is Gaussian and stationary. In this case, the noise back- 
ground for a single detector is Poisson distributed, with rate 
of events louder than pi is given by 



(22) 



where C depends upon the number of trials. With two detec- 
tors, the combined SNR is given by = pi + P2- Further- 
more, it is necessary to impose a single detector threshold on 
the SNR, denoted px- In this case, the Poisson rate of back- 
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FIG. 1: The distribution of the SNR of the loudest event. These dis- 
tributions were obtained by generating 100 time-shift loudest events 
from the distributions in Eq. l i24b with an event rate given by 
Eq. l l23b . with the normalization constant chosen so that the expected 
rate is unity at an SNR of 10. The upper plot shows the probability 
distribution while the lower plot shows the cumulative distri- 

bution Pb(p). The distribution is, as expected, peaked close to SNR 
of 10. The cumulative probability tends to zero at small SNR — it 
is almost guaranteed that there will be an event louder — and unity 
at large SNR — it is very unlikely to have an event louder than this. 
The dashed line in both plots shows the simulated result with an SNR 
of 9.95. 



ground events is 



p{p)^C{l + py2-pj,)e-p'-/' 



(23) 



The constant C can be determined from the expected loudest 
event. Given the Poisson rate p.{p), it is straightforward to 
calculate the distributions of Pb, pb and ub- 

Pb{p) = e-^^P'^ , pb{p) = \p'{p)\e-P^'^ andn^ = \p'{p)\ . 

(24) 

Given the background of ( l22b . at the expected loudest event, 

nB{Pm) ~ Pm- 

In our example, we choose pT = 5.5 and select C such that 
the expected loudest event has an SNR of 10, i.e. /i(10) = 1, 
equivalent to an SNR of about 7 in each detector For this 
search, we simulate 100 time shifts and obtain the loudest 
event for each. The distributions are plotted in figure [T] The 
features in these plots are due to the finite number of time- 
shifts performed, which lead to uncertainties in the recon- 
struction of the distributions. In addition, we simulate the out- 
put of the search, and obtain a loudest event with pm = 9.95 
which yields values of pb = 2.7, Pb = 0.25 and ub = 10.9. 



B. Calculating the Foreground 

The upper limit derived from a search depends upon the to- 
tal luminosity Cl to which a search is sensitive. This must 
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be evaluated at the SNR of the loudest event CL{pm)- In sec- 
tion |II] we have shown that the strength of the gravitational- 
wave signal at a detector depends primarily upon the effec- 
tive distance Dcff and chirp mass of the signal. In ad- 
dition, the measured SNR of signal with a given Deff and 
M. will fluctuate depending upon the noise in the detectors 
at the time. Thus, for each value of Z?off and Jv[, we can 
calculate the probability that the signal is observed with an 
SNR greater than p. This quantity is known as the efficiency, 
e(/9, Uoff : -^)- Since the sensitivity of a detector varies over 
the course of a search, the efficiency is measured by perform- 
ing Monte Carlo simulations. For a search involving several 
detectors, the efficiency will depend upon the effective dis- 
tance to the source from all detectors, which we denote in 
bold as Dcff. The efficiency calculation is discussed in Sec- 
tion HHbI] 

To calculate the total luminosity, we also require a measure 
of the blue light luminosity iB(I?cff, A^) as a function of Dcff 
and A^. This is calculated from the known luminosity density 
in the universe. For the initial LIGO detectors — with sensi- 
tivity to binary neutron star coalescences up to 40 Mpc — it is 
necessary to take into account the inhomogeneity of the local 
universe. This requires the construction of a catalog of nearby 
galaxies (see Ref. for details on how this is constructed). 
Armed with a galaxy catalog and a mass distribution for the 
binaries, the method of calculating Lb is described in Section 
1111 B 21 

Finally, given the efficiency and luminosity functions, the 
cumulative luminosity is given by 



Cl{p)^ J d-D,sdMe{p,-D,s.M)LB(Dcif,M). (25) 
Details of this calculation are given in Section lTllB 31 

/. Efficiency 

The efficiency of a search is evaluated by adding simulated 
inspiral signals to the data stream and evaluating the probabil- 
ity that signals with a given set of parameters are detected with 
SNR greater than pm- By performing a host of injections, it 
is possible to evaluate the efficiency as a function of both the 
chirp mass and effective distance. 

Figure 12] shows the simulated efficiency of the initial LIGO 
detectors, at an SNR of 7, as a function of the binary's effec- 
tive distance and chirp mass. The distance to which events 
are detectable increases with the chirp mass. This follows im- 
mediately from Eq. (|6]l for the inspiral waveform, from which 
we see that the amplitude is proportional to A4^^^ / Dcs- For 
binaries with a total mass less than 10 Mq, the inspiral stage 
of the evolution will sweep right across the sensitive band of 
the detector. Therefore, for low-mass signals, the detectability 
of a signal at a given detector will be dependent only on the 
chirp distance Dc of the binary: 




Dr := D. 



off 



/ Ml.4 



V M 



5/6 



(26) 
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FIG. 2: The efficiency of the initial LIGO detectors to coalescing 
binaries, as a function of the effective distance and chirp mass of 
the source. The efficiency curve is evaluated at an SNR of 7 in the 
detector. 



For higher mass signals, the coalescence will occur within the 
sensitive band of the detector, as can be seen from Eq. (|7]). 
In this case, there is no simple relationship between the effi- 
ciency, chirp mass and distance, and both mass and effective 
distance must be taken into account when calculating the effi- 
ciency. 

Let us now generalize to the case of several detectors. A 
given binary coalescence will have a different effective dis- 
tance in each detector, and therefore the efficiency will be a 
function of the effective distance in each detector In figure 
|3]we plot the efficiency against effective distance in the HI 
and LI detectors, for a 1.4-1.4 A/© binary. At small effec- 
tive distance in both detectors, the efficiency is unity, while 
at large distances it goes to zero. The shape of the constant 
efficiency contours depends upon two factors: the single in- 
strument threshold and the combined threshold. At large ef- 
fective distance in one instrument, but small in the other, the 
single instrument threshold limits the efficiency. For compa- 
rable distances, the efficiency is determined by the combined 
SNR. 



2. Astrophysical Model 

In order to calculate the total luminosity to which a search 
is sensitive, we require the luminosity density as a function 
of the effective distance and mass. This involves combining 
a model of the location and luminosity of galaxies with the 
antenna response functions of the detectors, given in Eq. (O. 

Let us consider a source from a given galaxy. The effective 
distance to the source depends upon the physical distance, and 
four sky angles — the sky location relative to the detector, the 
inclination and polarization angles. Equivalently, the location 
of the source can be described by A = {D,a,d, L,ip,t), where 
the right ascension a, declination S and sidereal time t serve 




10 20 30 40 50 60 
DeHLIIO(Mpc) 



FIG. 3: The detection efficiency as a function of the effective distance 
for two detectors. The efficiency at small distances is unity, while at 
large distances it is zero. The transition is governed by two effects. 
At large distances in one detector, the efficiency is limited by the sin- 
gle detector threshold of 5.5. At large distances in both instruments, 
the efficiency is determined by the combined SNR threshold of 10. 
The variations in the contours are due to the fact that a finite number 
of injections — 100 in each bin — were used when generating the 
efficiency curve. 
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LIIO(Mpc) 



FIG. 4: The luminosity distribution in the nearby universe as a func- 
tion of the effective distance in the Hanford and Livingston detectors. 
The distribution is obtained by assuming that the binary coalescences 
in a given galaxy are dependent upon the blue light luminosity of that 
galaxy and are uniformly distributed in sidereal time, inclination and 
polarization. The color bar indicates the available blue light lumi- 
nosity (Z/io/Mpc^). The number increases with distance in both 
detectors and is greatest on the diagonal. The off diagonal spread 
in luminosity is due to the different alignments of the Hanford and 
Livingston detectors. 



to specify the source sky location. For sources from a par- 
ticular galaxy, the distance, right ascension and declination of 
the galaxy {Di, ai, Si) are known. Then, assuming that binary 
coalescences are uniformly distributed over the sidereal day, 
inclination and polarization angles, we obtain a distribution 
for the observed effective distances of sources from a given 
galaxy 



PiiDcs) 



dXp{t)p{i) p{ip) 5{D ~ A) 5{a - a^) 



(27) 



where = 1/ (sidereal day), ^(t/') ~ l/27r and p(t) = 
sin(i)/2. In practice, this distribution is most easily obtained 
by simulating many signals, at random times and orientations, 
from a given galaxy. 

In Ref. Isl], the compact binary coalescence galaxy (CBCG) 
catalog has been compiled to a distance of 100 Mpc. For 
galaxies in this catalog the sky location and distance to the 
galaxy are known. In addition, the apparent magnitude niB.i 
in blue light of the galaxy is measured. The luminosity of the 
galaxy is obtained from its distance Di and apparent magni- 
tude as 



Ls.i ( Di \' 



L 



(28) 



where Lb,q = 2.16 x lO^^ergs/sis the blue solar luminosity, 
and Mb,q is the (absolute) blue solar magnitude j34]. 

Given the distribution of effective distances for each galaxy, 
it is straightforward to obtain the total luminosity as a function 



of effective distance by summing over all galaxies: 



(29) 



As before, this can be easily generalized this to a distribution 
of luminosity as a function of effective distance for several 
detectors. In Figure |4] we make use of the galaxy catalog of 
lit] to plot the luminosity as a function of the effective distance 
at the Hanford and Livingston detectors. 

Finally, we must include the mass distribution p(A^). The 
available luminosity is expressed as a function of effective dis- 
tance and chirp mass as 



LBiD,ff,M) = (^LBaP^{D,^f?j 



piM) . (30) 



For neutron star binaries, the mass distribution can be taken 
from observations Issll or, alternatively, from population syn- 
thesis i^. For binaries containing at least one black hole, the 
lack of observations leads to greater uncertainty in the mass 
distribution, whence it is more natural to calculate a mass de- 
pendent rate limit. 



3. Calculating the cumulative luminosity 

The cumulative luminosity available to a search is obtained 
by multiplying the luminosity distribution in Eq. (|30] l against 
the efficiency e(p, Deft, A^) and integrating over the mass and 
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effective distance: 



Cl{p) = J ciDeff dM e{p, Deff , M) Lb(D,s) p{M) . 

(31) 

For a given mass, this reduces to multiplying the data from 
Figures [3] and |4] and integrating over the effective distances. 
This must then be integrated over of the chirp mass to obtain 
the cumulative luminosity. For low-mass systems, the cal- 
culation is greatly simplified by recalling that the efficiency 
depends only upon the chirp distance (|26] | at a given site. 

In a similar manner, the derivative of the cumulative lumi- 
nosity C^((0„j) can be calculated. This is done by calculating 
the rate of change of efficiency as a function of SNR, evalu- 
ated at With this information, we can calculate: 



np = 



\C'L{Pm)\ 
CL{pm) 



(32) 



In cases where the mass distribution is not known, the lu- 
minosity can be expressed as a function of the mass: 

Cl(p, M) = j dDeff e(p, Deff, M) Lb{T>,s) . (33) 

Then, the rate limit is naturally expressed as a function of the 
mass. 

We can compute the for the example with a loudest event 
SNR of 10. To simplify the calculation, assume that the mass 
of all binary components is IAMq. Then, the total luminosity 
and its derivative are 



Cl(p™) = 540Lio and Ci(p,„) = 120L 



10 ■ 



(34) 



This gives a value of np = 0.22. For comparison, we can 
calculate the expected value for a uniformly distributed popu- 
lation. In this case, Cl oc p~^, whence jij? w 3/p„i, which is 
similar to the calculated value. The difference is due to the fact 
that at the distances under consideration, the non-uniformity 
of the luminosity distribution is important. 



C. Obtaining an Upper Limit 

In the preceding sections, we have described all of the 
pieces which are required in calculating an upper limit from a 
search for gravitational waves from binary coalescence. The 
calculation of the cumulative luminosity Cl depends upon the 
efficiency of the search, and the astrophysical distribution of 
signals. We have argued that in the non-spinning case, these 
are both described as functions of effective distance and chirp 
mass. The likelihood A depends upon the foreground and 
background distributions, encoded in np and riB- The first 
of these, n p is obtained from the cumulative luminosity and 
its derivative, while ub is, determined from analysis of events 
arising in analysis of time shifts of the data. 

The formula for the upper limit was given in Eq. ( |2TI ). Tak- 
ing the prior distribution p{R) of the rate to be uniform, we 
obtain a Bayesian upper limit with confidence a as 



A 



1 + A 



RTCl {Pn 



(35) 



For a representative loudest event, we obtain ns ~ Pm and 
np w i/ Pm- Therefore, A ~ "^1 P^n- For ™y realistic 
search, with a loudest event of SNR around 10, this implies 
that A ^ 1, and therefore the loudest event is most likely 
background. To obtain a mass dependent upper limit, this 
analysis is repeated for different values of 7W making use of 
the mass dependent luminosity function CL{pm, -M) to obtain 
a rate limit R{Ai). This method neglects the fact that the 
distribution of background events is also mass dependent by 
using the same loudest event for all masses. 

Plugging in the values obtained in the previous sections, 
we have A = 0.22/10.9 = 0.020, and Cl(/9,„) = 540 Lio. 
Assuming a year of analysis time, the limit would be 



R. 



2.35 



90% 



CpiPm) T 



= 4.3 X IQ-^L-iyr-i 



(36) 



This gives a reasonable estimate of the expected BNS upper 
limit in the absence of a detection in the LIGO S5 science run. 

How does this compare with astrophysical predictions? The 
rate of galactic binary neutron star mergers is 8.30l;g*'g^^ x 
10~^?/r^^, at 95% confidence ll36ll . Assuming the rate is in- 
deed a function of the blue light luminosity alone, and using 
1.7^10 as the blue light luminosity of the Milky Way |f36|l, 
gives an optimistic rate of 1.7 x 10~'^L^q yr^^, which is a 
factor of 25 lower than the expected upper limit from the anal- 
ysis of one year of data from initial LIGO. 



IV. SYSTEMATIC ERRORS 

In the previous section, we have described a method for 
calculating an astrophysical upper limit or interval for the rate 
of binary coalescences from the results of a gravitational-wave 
search. The probability distribution for the rate is dependent 
on four quantities: the prior distribution p{R), the cumulative 
luminosity Cl, the likelihood A of the loudest event being a 
signal, and the analysis time T. Of these quantities, only the 
analysis time can be unambiguously measured. The choice of 
prior distribution p{R) will affect the posterior distribution. 
However, we take the prior as an input to the analysis and 
do not consider uncertainties associated to the choice of prior 
There are numerous systematic errors which will affect the 
measured values of both the luminosity and likelihood. These 
systematic effects can be broadly split into five categories: 

• Imprecise knowledge of the astrophysical distributions 
of the mass and distance of binaries. 

• Differences between the physical signal and the non- 
spinning, restricted post-Newtonian waveforms. 

• Statistical fluctuations in the measured efficiency. 

• Uncertainties in instrumental calibration. 
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• Errors in the calculated likelihood A, arising from the 
above uncertainties and errors in the background esti- 
mation. 

In this section, we describe the various sources of uncer- 
tainty and analyze their effect. Additionally, we perform a 
marginalization over these uncertainties to produce the rate 
distribution. 



A. Uncertainties in population model 

The cumulative luminosity of a search will depend criti- 
cally upon the astrophysical model used. In particular, the 
luminosity distribution is sensitive to both the location and 
luminosity of galaxies within the reach of the search. In addi- 
tion, since the amplitude and frequency range of gravitational 
waves emitted by a binary coalescence is mass dependent, the 
cumulative luminosity will also be dependent upon the astro- 
physical mass distribution. In Section ITV A II we discuss the 
systematics associated to uncertainties in galaxy distribution, 
while in llV A2l we investigate the effect of changing the mass 
population. 



1. Galaxy Distribution 

The sky position of nearby galaxies is known accurately 
enough that errors in the right ascension and declination will 
not affect the cumulative luminosity. However, the distances 
and luminosity of galaxies are not so well known, whence 
these uncertainties must be taken into account when calcu- 
lating the total luminosity. Indeed, the luminosity of a galaxy 
is not directly measurable, instead it is calculated from the ap- 
parent blue magnitude niB.i and distance Di using ( l28T l. The 
uncertainties in distances vary from less than 10%, from ob- 
servations of Cepheids in nearby galaxies observed with the 
Hubble Space Telescope, to uncertainties up to 30% for more 
distant galaxies. Additionally, there are uncertainties in the 
apparent magnitude ofgalaxies which vary from AmB,i = 
0.3to AjTiB,, 0.38 il. 

We begin by considering an error in distance of ADi to a 
galaxy at distance Di. The change in distance of the galaxy 
will lead to a linear change in the effective distance of all 
sources from that galaxy. More precisely, changing the dis- 
tance from Di to {Di + ADi) will change to distribution of 
effective distances fmmpi{Dcs) to pA,j:(^ofr) where 

PA,,:[^cff(l + AD,/D,)] = p,{D,s) . (37) 

Thus, the average effective distance to a source will increase 
by (1 + ADi/Di). This will not have any effect for close 
(or very distant) galaxies where the efficiency to sources from 
that galaxy is unity (or zero). However, for galaxies at dis- 
tances where the efficiency is rapidly changing, this can be a 
significant effect, reducing the efficiency when the distance to 
the galaxy is increased. 



Since the luminosity of a galaxy is inferred from its mea- 
sured distance and apparent magnitude, a change in the dis- 
tance will also affect the calculated luminosity. It follows di- 
rectly from ( l28T l that a change in distance of ADi, leaving the 
apparent magnitude unchanged, will yield a change in lumi- 
nosity ALB,i given by 

Lb,, + ALb^, _ f D., + AD, V ^^g^ 
Lb,i V A / 

Thus, the inferred luminosity will increase if the distance to 
the galaxy increases. This effect will be significant for all 
galaxies to which the search has non-zero efficiency. 

It is also straightforward to calculate the effect of errors in 
the reported apparent magnitude. We have from Eq. ( |28T l 

LB,i + ALb^i ^ -^Q(AmB,,/2.5) _ (39^ 

Lb,! 

Therefore, an increase (decrease) in magnitude will lead to an 
increase (decrease) in Cl. The reported errors in the CBCG 
catalog are Am^ ^ between 0.3 and 0.38, giving correspond- 
ing uncertainties in the luminosity of 32% and 42% respec- 
tively. 

In order to estimate the uncertainty in the luminosity Cl, it 
is necessary to vary the distance and magnitudes of all galax- 
ies for which the search has non-zero sensitivity. We expect 
that the errors are independent for different galaxies, although 
it is certainly conceivable that there is an overall systematic 
for a given galaxy catalog. However, to be conservative, we 
calculate the errors by moving all galaxies either closer or fur- 
ther As with the distance error, we take the conservative error 
by increasing or decreasing the magnitude of all galaxies to- 
gether 

At large distances, the luminosity distribution can be taken 
as homogeneous and isotropic. In this case, the blue lu- 
minosity density is calculated directly from observations as 
Lb = (1.98 ±0.16) x lO-^Lio/Mpc Thus, at distances 
greater than 40 Mpc, the uncertainty in luminosity can be cal- 
culated directly. 

Recent papers have suggested that due to the large coales- 
cence times for binary inspiral, a significant fraction of them 
might occur in elliptical galaxies where the star formation 
rate is low. The general loudest event formalism presented 
in Ref. jsj can be used in this case, and one might envision 
introducing two unknown rate parameters Rb which depends 
upon the blue light luminosity discussed above and Re which 
is a second contribution due to elliptical galaxies. Then, the 
rate would depend upon two parameters, and a given search 
could be used to place a confidence bound in the two dimen- 
sional space spanned by them. If the corrections from includ- 
ing elliptical galaxies prove to be significant, this effect will 
be folded into future rate calculations. 



2. Binary mass distribution 

There is significant uncertainty in the mass distribution of 
coalescing binaries. Several binary neutron star systems, and 
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1.0 L2 1.4 1.6 1.8 2.0 2.2 

Chirp Mass (Mr,) 




Effective Distance (Mpc) 

FIG. 5: Effect of the mass distribution on the search efficiency. The 
upper plot shows two different binary neutron star mass populations: 
one taken directly from observed Neutron star masses, which are fit- 
ted to a Gaussian mass with a peak at 1.35A/0 and width of O.O4M0 ; 
the other from population synthesis models in Ref. The lower 
plot shows the efficiency as a function of effective distance for the 
initial LIGO detector at SNR 7, where mass has been folded in. 

significant numbers of single neutron stars, have been ob- 
served as pulsars, allowing us to place some restrictions on the 
mass distribution. However, the mass distributions presented 
in (for example) |@| are produced using large scale simula- 
tions which must assume an equation of state for the nuclear 
material. For stellar mass black hole binaries, there is little 
restriction on the mass distribution. These uncertainties lead 
us to place mass dependent upper limits. However, it is still 
instructive to look at the sensitivity of our binary neutron star 
search for various mass distributions. 

As discussed above, the distance to which a source can be 
observed is dependent upon its chirp mass. To leading or- 
der, the amplitude of the emitted gravitational radiation, and 
hence the distance to which the source can be observed, is 
proportional to A4^^^. The astrophysical mass distribution of 
neutron stars can have a significant effect upon the distance to 
which sources are observable. For a 1.4 — 1.41/0 solar mass 
binary (M ~ 1.22A/q), the 50% efficiency point for initial 
LIGO at SNR 7 occurs at 40 Mpc, whereas for a 3.0 - S.OA/q 
binary that is increased to 75 Mpc. So, an astrophysical popu- 
lation containing higher mass binaries will increase our range. 
As an example, let us consider two mass distributions; 

• The distribution of observed masses of radio pul- 
sars issll . namely a Gaussian distribution peaked at 
1.35M0, with a width of 0.04Mq. 

• The distribution from Ref. obtained from population 
synthesis models. This is the mass distribution which 
was used in obtaining the LIGO S 1 and S2 results given 
inRefs. SlU. 

Figure |5] shows the distribution of chirp mass for the two dif- 
ferent cases. In both cases, the peak is at around the same 



chirp mass of 1.2Mq, corresponding to binaries with com- 
ponent mass around 1.4M0. However, there is a significant 
fraction of higher mass neutron stars in the population synthe- 
sis distribution. 

Given a mass distribution, we can integrate the efficiency 
over the mass to obtain efficiency as a function of effective 
distance alone, 

e{D,s)^ j dMp{M)e{D,s,M). (40) 

This curve is also plotted in Figure |5] for the two mass distri- 
butions of interest. Varying the mass distribution has a 20% 
effect on the distance at which the efficiency reaches 50%. 
Since the choice of mass distribution can have a significant 
effect on the efficiency, in ifioll the upper limits from LIGO 
searches are reported for a stated sample distribution. The 
mass uncertainty is not folded into the systematic errors on 
the upper limit. 

B. Uncertainties in the Waveform 

In several places, we have assumed that both the phys- 
ical waveforms and the search templates are described by 
the post-Newtonian approximation and we have ignored the 
effects of spin. The methods presented here are not tied 
to the particular templates used and extend immediately to 
searches using other types of waveform, or even to untem- 
plated, excess-power type analyses iT2tl . The only require- 
ment is that it is possible to meaningfully associate an SNR to 
both noise and simulated signals, and hence obtain a loud- 
est event. Throughout, we have also modeled the true as- 
trophysical waveforms as those obtained from the restricted, 
non-spinning, post-Newtonian calculation. This can have a 
significant effect on the interpretation of the result. 

In order to explore how the uncertainty in the simulated 
waveforms will affect our result, it is useful to introduce the 
notion of match. For a given set of parameters (in particu- 
lar the component masses and effective distance), we denote 
the true waveform by ht and those used in simulations as hg. 
The difference between the true and simulated waveforms is 
encoded in the match M defined as 

M = . (41) 

^{ht\ht){hs\hs) 

If the waveforms agree perfectly, the match will be unity, 
while in all other cases it will be less than one. 

Differences between the post-Newtonian approximation 
and the true waveforms have been examined in some detail. 
In Refs. ll39il40ll . this has been done by comparing the post- 
Newtonian waveforms to those obtained from black hole per- 
turbation theory. The results indicate approximately a 10% 
loss of SNR (i.e. a match of 90%) due to inaccurate modelling 
of the waveforms for low-mass systems, with the effect be- 
coming more pronounced for higher mass ratios. In Ref. ll4lll 
a similar result was obtained by comparing waveforms at dif- 
ferent post-Newtonian order Recent breakthroughs in numer- 
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ical simulations of black hole and neutron star coalescences 
will allow for further exploration of this issue. 

In order to fully address the issue of spin, it will be neces- 
sary to perform substantial Monte-Carlo simulations of spin- 
ning waveforms. For the time being, we rely on estimates 
provided by Apostolatos [[Till in which he shows that less than 
10% of all spin-orientations and parameters consistent with 
binary neutron stars provide a loss of SNR greater than 5%. 

In BNS searches, where the same waveforms are used as 
templates and simulated signals, it is straightforward to cal- 
culate the effect of any mismatch. In this case, the mismatch 
between the astrophysical waveform and the post-Newtonian 
approximation will lead to an over-estimation of the observed 
SNR from a given binary coalescence. Specifically, if pt and 
Ps are the SNRs associated to the true signal and simulation 
respectively, then 



Pt = PsM^ 



' {ht\ht) 
{hs\hs) 



(42) 



Therefore, the SNR associated to a true signal is reduced by a 
factor M from what is observed in a simulation. 

There is an important difference between the waveform un- 
certainties and the other systematic errors discussed in this 
section. In the simulation, we are using the same waveform 
for injection and detection. In reality, the true astrophysical 
waveforms will not match precisely the detection family. This 
will lead to a decrease in the overlap between the astrophysi- 
cal and detection families. It is not possible for this to lead to 
an increase as the match cannot be greater than unity. Thus, 
the waveform errors can only serve to decrease the cumulative 
luminosity available to a search. However, in cases where the 
simulated signals and templates do not match precisely, it is 
possible that the "true" waveforms will have a better match 
with the templates than the simulations do. Then, errors in the 
waveform may cause us to underestimate the efficiency of the 
search. 

Returning to Eq. ( l42b . we see that the overall normaliza- 
tion of the waveform will also affect the SNR. Generally, it is 
assumed that the amplitude of the simulated and astrophysical 
waveforms are in good agreement, namely {hs\hs) ~ {ht\ht). 
In Refs. ll42ll43ll . the authors consider the effects of higher or- 
der post-Newtonian corrections to the amplitude. These are 
particularly important for higher masses, especially when the 
ratio of the component masses is large. Furthermore, in ll42ll 
it has been noted that the inclusion of amplitude corrections 
actually decreases the overall amplitude of the signal. Thus, 
even though neglecting amplitude corrections may not signif- 
icantly affect the detection ability of a search, it can still have 
an effect on the interpretation of results. This effect is not im- 
portant for BNS systems, but does become important in higher 
mass, asymmetric systems. 



C. Uncertainties in the instrumental calibration 



a simulated signal differs from the SNR pt that would be as- 
sociated with a real signal with the same parameters due to 
uncertainties in the instrumental calibration. 

In calculating the efficiency of a search, simulated events 
are added in software to the data after it has been recorded. 
Therefore, any uncertainty in instrumental calibration will not 
affect the software injections in the way it will a real signal. 
To quantify this effect, we focus on the differences between 
a true signal with given masses and effective distance _Deff 
and a simulated signal with the same parameters. To simplify 
matters, we assume that the waveform exactly matches one 
of the search templates (i.e. ignore the waveform uncertainty 
discussed above), in which case 



(43) 



where hdf) is the waveform introduced in ([Sjl, (j)o is an 
arbitrary phase, and n{f) is the detector noise. Then, by 
substitution into Eq. (fT2l i. it is straightforward to verify that 
= + 2 as expected. 
The output of a gravitational-wave detector is not the grav- 
itational wave strain, it is an uncalibrated signal v{t). This 
output is related to the strain by 



- Rtimf) 



(44) 



where Rt{,f) is the true response function of the instrument. 
The process of calibrating the data involves the construction 
of a response function R{f). Due to calibration uncertainties, 
the reconstructed response will differ from the true response. 
These calibration uncertainties can, to some degree, be inde- 
pendently tested by performing "hardware injections" into the 
detectors ll44ll . To calculate the effect of calibration uncertain- 
ties, let us follow 1I45I1 and write the measured response as 



R{f) = r{f)e 



(45) 



and the true response as 

Rt{f)^[r{f) + 6r{f)V^'^^'^''''^^'^^- (46) 

where 5r and 56 are the amplitude and phase parts of the cal- 
ibration error which are assumed to be small. 

In the event that a gravitational wave with strain h{t) im- 
pinges on the detector, the calibrated strain reconstructed from 
the output of the detector will be given by:^ 



W) = 



RtU) 



(47) 



The SNR for an event observed in the detector can be calcu- 
lated by substituting Eq. ( |47] i into the expression fT% for the 
SNR. Then, making use of the expressions ( |45] ) and ( |46l ) for 



When performing simulations of the gravitational-wave 
signal from a coalescing binary, the SNR ps associated with 



^ For simplicity, we ignore contributions from the noise in the following. 
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the response functions, and expanding in powers of br jr and 
59 v^Q obtain 



r? - f? 

Pt - Ph ^2 




df 



SnU) 



(48) 



Thus, the error is linear mSr but quadratic in 59. Furthermore, 
it follows directly from the Schwarz inequality that 



|(M6'|/ic)P < a'^{hc5e\K5e). 



(49) 



This guarantees that the magnitude of the second (positive) 
contribution to the SNR from the phase error is always less 
than the first (negative) one. Hence, an error in the phase of 
the response function will always serve to reduce the recov- 
ered SNR. 



1. Application to the LIGO instruments 

The calibration report from the LIGO S4 science run ||4C 
provides an in depth analysis of the calibration of the LIGO 
instruments. Here, we briefly review those details which are 
relevant for coalescing binary searches. 

In the LIGO instruments, the response function is deter- 
mined from the sensing function C{f) and the open loop gain 
function G{f) as 



Rif) 



l + Gjf) 
C{f) 



(50) 



These calibration functions, C(/) and G(/) are measured at 
intervals during an analysis. At intervening times, it is as- 
sumed that their functional form is unchanged. However, due 
to the changes in light power stored in the arms, a time depen- 
dent rescaling ^{t) of the both the sensing function and open 
loop gain is required. Thus, the response function for any time 
can be expressed as: 



R{f,t) = 



l + lit)Goif) 

l{t)Co{f) 



(51) 



The uncertainties in the various components of the response 
function are calculated in detail in Ref. ll46ll . These can be 
summarized by the statement that the uncertainties in the re- 
sponse function are of order 5% in the amplitude of the re- 
sponse function and 5° in phase. Substituting this into our 
general expression (l48T l, we obtain: 



Pt 



It 



\5r\ \5e\^ 



p [It 0.05]. (52) 



Note that the error in the measured SNR is primarily due to 
the amplitude calibration uncertainty. 



D. Uncertainties in the measured efficiency 

The effects of discreteness of the template placement, er- 
rors in the estimates of the power spectral density S{f) used 
in the matched filter in Eq. ( fT2] i. and trends in the instrumen- 
tal noise are all accounted for by the Monte-Carlo simulation. 
However, in the Monte Carlo simulation, only a finite num- 
ber of injections are performed as these are computationally 
costly. Thus, the efficiency plots, such as that shown in Figure 
12] have an associated statistical error. The efficiency plots are 
produced by binning up the parameter space and calculating 
the efficiency in each bin as 



Nf 



Nf + N„ 



(53) 



where Nf and Nm are the number of found and missed in- 
jections respectively. Then, assuming binomial errors for the 
efficiency, the variance in the efficiency is: 



NfN„ 



{Nf + N„ 



{Nf + N,n 



(54) 



As expected, the variance is inversely proportional to the 
number of injections performed. Furthermore, it is clear from 
(l54l i that the statistical uncertainty in the efficiency is greatest 
when the efficiency is close to 50%. The Monte Carlo error 
in the luminosity is obtained by multiplying the error in the 
efficiency by the luminosity:^ 

{ACL{p)f= <^'!{p,'DcS,M)Lb{B,s)p{M). 

(55) 

When performing a search on real data, software injections 
are computationally costly. From Eq. ( |55] l, we see that simula- 
tions are most efficiently performed when concentrated where 
the efficiency is close to 50% and there is a significant contri- 
bution to the luminosity. Furthermore, for low-mass signals, 
making use of the chirp distance (|26] l reduces the dimension 
of (ISST i by one and consequently reduces the size of the asso- 
ciated Monte Carlo errors. 



E. Marginalization over Uncertainties to obtain an Upper 
Limit 



In section |III] we have described how to calculate an up- 
per limit from an inspiral search, making use of the loudest 
event statistic. This involves calculating three quantities, the 
amount of time searched T, the total luminosity Ci,(p„i) to 
which the search is sensitive, and the likelihood ratio A of the 
event being foreground. In this section, we have discussed 



^ Strictly, there is a second Monte Cai'lo uncertainty due to errors in tlie cal- 
culated luminosity. However, the luminosity density of Figure |4] can be 
calculated to good accuracy with minimal computational cost so that these 
errors will be insignificant. 
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various systematic errors which affect our ability to measure 
these quantities. Here, we will describe how these uncertain- 
ties can be marginalized over to produce an upper limit which 
takes them into account. 

Let us begin by considering a single uncertainty (for ex- 
ample the distance error) which will effect the luminosity Cl- 
Then, by evaluating errors associated to this uncertainty, we 
obtain a probability distribution for the cumulative luminosity, 
PdiCh)- In order to marginalize over the distance uncertainty, 
we simply evaluate 

T,B) = j dCL PdiCL) p{R\Pm, T, B,Cl), (56) 

where the probability distributions on the right hand side must 
be normalized to unity. This is straightforward when there is 
only one error to take into account. However, in the preceding 
sections, we have detailed several errors. It is not practical 
to integrate over all of these errors independently, so we per- 
form a Monte Carlo integration. For the majority of the errors, 
we use a Gaussian with standard deviation given by the value 
stated above (and truncated so that the total Lio will never 
be negative). Since the magnitude error affects the luminos- 
ity exponentially, assuming a Gaussian error on this leads to 
a log-normal distribution for the luminosity. Finally, for the 
waveform error, we use a 1 -sided Gaussian, i.e. one which 
can only decrease the cumulative luminosity. Then, the Monte 
Carlo integral is performed by sampling many times from the 
appropriate distributions. 



F. Uncertainties in tlie likelihood 

We have cataloged various uncertainties which will affect 
the calculated luminosity Cl for a given search and shown 
how they can be marginalized over to obtain a final distribu- 
tion for the rate. In addition, we need to examine the effect 
any uncertainty in the estimation of the likelihood A, will have 
upon the upper limit. To do so, we will once again marginalize 
over this nuisance parameter (as in Eq. (|56]l), to obtain a distri- 
bution which is independent of A. However, due to the simple 
manner in which A enters the probability distribution ([TtT i for 
the rate, to leading order the marginalization procedure has 
no effect on the distribution More precisely, provided the 
uncertainties in A are small compared to (1+A), the marginal- 
ization procedure serves to replace A by its expectation value. 
We have argued that for a typical search where the loudest 
event is consistent with the background, we expect A <C 1 
whence it is safe to neglect uncertainties in its value. 

There are cases when the estimated background from time 
shifts will not accurately reflect the background. An obvi- 
ous example of this is when computing the background for 
the co-located Hanford detectors. It is well known that there 
are correlated noise sources which will produce inspiral trig- 
gers in the two instruments simultaneously. While every ef- 
fort is made to remove these times by examination of aux- 
iliary channels, and in particular the seismic data, there are 
still invariably some correlated noise event which occur in the 



Hanford data. In this case, time shifts will not give a good es- 
timate of the background. There are other ways to estimate the 
background, such as using "reverse chirp" filters (obtained by 
time inverting the template). However, in calculating an up- 
per limit, an underestimation of the background will lead to a 
conservative upper limit being placed. 



G. An Example 

To illustrate the issues associated with these systematic un- 
certainties, we will calculate them for the example introduced 
earlier For our example, the simulated loudest event had a 
combined SNR of pm = 9.95, which corresponds to a single 
instrument SNR around seven. At this SNR, the 50% effi- 
ciency for BNS occurs at 40 Mpc. The cumulative luminosity 
of such a search is CL{pm) = 540 iio. Finally, the value of 
the likelihood for our simulated results is A = 0.02. Now, 
we turn our attention to the systematic uncertainties discussed 
above. We will evaluate the effect of each of these on the 
cumulative luminosity. 

Distance uncertainties are obtained by moving all galaxies 
either closer or further away, by the appropriate fraction given 
in the CBCG-catalog. This changes the effective distance to 
sources according to Eq. ( [37] i and the luminosity according 
to Eq. ( l38T l. This yields a change in cumulative luminosity 
of 90iio with a luminosity decrease as galaxy distances are 
increased. Uncertainties in the magnitudes of galaxies are 
taken into account by rescaling their luminosities according 
to Eq. (|39] l and lead to an error in the cumulative luminosity 
oflOOiio. 

Based on the discussion of Section |TV B| for BNS systems 
we choose to use a 10% uncertainty in the astrophysical wave- 
form. This is simulated by reducing the efficiency of both de- 
tectors by the amount. In other words, we rescale the axes in 
Figure [3] downwards by 10%. The effect in our simulation is 
a change in the luminosity of 160 Liq. 

The calibration uncertainty is calculated by varying the ef- 
ficiency curve accordingly, in a similar manner to that used 
for the waveform errors. Since the calibration errors in dif- 
ferent instruments are independent, we consider them one at 
a time. In our example, we take a 5% calibration uncertainty 
in both the HI and LI detectors. This is calculated to lead to 
a 37 LiQ variation in the luminosity in HI and a 34 Liq varia- 
tion in LI. These numbers are very similar, as expected. They 
differ slightly due to the fact that certain galaxies are better 
aligned for one detector than the other 

The Monte Carlo uncertainty can be calculated according 
to Eq. (|55T l. As an example, Figure|2]was produced with 100 
injections in each bin. This gives a Monte Carlo error of 3 iio 
which is already weU below the errors due to astrophysical and 
waveform uncertainties. 

The probability distribution for the cumulative luminosity 
taking into account these systematic errors is shown in Fig- 
ure|6] To generate the distribution, the uncertainties described 
above are used as the standard deviation of Gaussian distribu- 
tions sampled via a Monte Carlo process to obtain the result. 
As is clear from the figure, the width of the probability dis- 
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FIG. 6: The distribution of the cumulative luminosity after marginal- 
izing over the systematic errors. The histogram was produced by 
using 100,000 samples, and includes all the noise sources described 
in this section. The dashed vertical line gives the luminosity before 
marginalization. The fact that the peak of the distribution is below 
this is due to the fact that the waveform errors are one sided. 
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FIG. 7: The cumulative distribution on the rate being greater than a 
given value. We are interested in the 90 % upper limit and therefore 
take the rate for which the distribution is equal to 0.1. The distri- 
bution is plotted for both the original luminosity of 540Lio and the 
marginalized luminosity distribution shown in figure |6] The large 
tail on the marginalized rate is due to the width of the luminosity 
distribution. 



tribution is significant. The greatest contributions to the un- 
certainty come from the astrophysics — the uncertainty in the 
luminosity distribution of nearby galaxies — and uncertainties 
in the waveform. The systematic effects due to instrumental 
calibration and Monte Carlo errors are sub-dominant. 

It is illustrative to consider the magnitude error and the 
waveform error in more detail. Although these errors are sim- 
ilar in magnitude, the waveform error is one sided, and thus 
has a more significant effect on the result. Although the mag- 
nitude error is 100 Lio, marginalizing over it only increases 
the 90% confidence upper limit from 4.3 x L^^^yr^^ 
to 4.5 X IQ^'^ L^^yr^^ . In contrast, marginalizing over 
only the waveform error increases the upper limit to 6.2 x 
10~^L^^yr~^ . The effect is much more significant, even 
though the magnitude of the two systematic errors is similar 
The reason is that the waveform error can only decrease the 
sensitivity of the search. Indeed, by modeling the waveform 
error as a 1 -sided Gaussian, we are significantly changing the 
mean value of CL{pm), reducing it from 540 Lio to 410 iio, 
which accounts for a large fraction of the increase in the re- 
ported upper limit. When we include all of the systematics 
together, we obtain a luminosity distribution as shown in fig- 
ure |6] 

Finally, we can make use of the above distribution for the 
cumulative luminosity in order to construct the posterior prob- 
ability distribution for the rate. We do this both for the un- 
marginalized and marginalized cases, and these are shown in 
Figure|2] This shows the cumulative posterior probability dis- 
tribution for the rate of binary coalescence given the search. 
The figure shows both the un-marginalized and marginalized 
distributions. The marginalized distribution takes into ac- 
count the systematic errors discussed above and falls off more 
slowly due to the width of the luminosity distribution from 
Figure |6] The final, marginalized upper limit for this example 



can be read off from Figure|2]as 7.3 x 10 ^ L^^ yr ^. 



V. SUMMARY AND CONCLUSIONS 

The astrophysical interpretation of the results is a critical 
part of any search for gravitational waves. In this paper, we 
have described a method to obtain an astrophysical rate up- 
per limit or interval from the results of a search for coalescing 
binaries. This method has been used in obtaining the results 
from the LIGO S3 and S4 science runs llOll . An astrophys- 
ical interpretation of a search result requires a good under- 
standing of both the detector sensitivity and the relevant astro- 
physics. We have argued that the detector sensitivity can be 
expressed in terms of the efficiency and furthermore that, for 
non-spinning systems, this efficiency is dependent only on the 
chirp mass and effective distance of the binary relative to the 
detectors. Furthermore, for low-mass systems, these can be 
combined into a single quantity, the chirp distance ( |26] |, which 
characterizes the amplitude. The use of effective distance is 
only appropriate for non-spinning binaries. Once spin is in- 
cluded, the orbital plane of the binary precesses whence the 
effective distance is not a constant over the course of inspiral. 
Despite this, the methods presented here could be extended to 
spinning binaries by making use of the "expected SNR" {ph 
in Eq. ( fT4l i) to characterize the amplitude of the signal. How- 
ever, while it is clear that the efficiency for a non-spinning bi- 
nary will be a function only of the effective distance and chirp 
mass, it is not obvious that the expected SNR and masses will 
completely characterize the efficiency to spinning systems. 

The relevant astrophysics is encoded in the expected distri- 
bution of coalescing binaries in the universe. Following |@], 
we make the assumption that compact binaries are distributed 
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according to the blue light luminosity. Making use of a cat- 
alog of nearby galaxies (such as the one in Ref. ISfl), we ob- 
tained an expression for the total luminosity to which a given 
search is sensitive. The loudest event statistic allows us to 
obtain a probability distribution for the rate of binary coales- 
cence given the results of a search. This distribution depends 
upon the cumulative luminosity discussed above as well as an 
understanding of the rate of background, noise events present 
in the data. Finally, the posterior distribution for the rate can 
be used to calculate an upper limit or rate interval for the oc- 
currence of binary coalescence. 

There are numerous systematic uncertainties involved in the 
calculation of the rate which we have discussed in detail. The 
dominant errors arise due to uncertainties in the distribution of 
nearby galaxies and imprecise knowledge of the gravitational 
waveform emitted by coalescing binaries. It is reassuring that 
the errors associated with our understanding of the detectors 
and the analysis, such as calibration and Monte Carlo statis- 
tical uncertainties, are less significant than the physical and 
astrophysical uncertainties discussed above. 

Uncertainties in the measurement of distances and appar- 
ent magnitudes of nearby galaxies leads to an uncertainty in 
the total luminosity available to a search. This in turn af- 
fects the reported rate limit. Additionally, the unknown mass 
distribution of coalescing binaries will have a significant ef- 
fect on the reported rate. To mitigate this effect, rates are re- 
ported as a function of mass. In the future, gravitational-wave 
observations of binaries by advanced detectors will provide 



improved knowledge of both the mass and location distribu- 
tion of binaries. Although the post-Newtonian waveform is 
known to a high level of accuracy, uncertainties in the wave- 
form, the neglect of spin and amplitude corrections lead to 
significant uncertainties in the sensitivity of the search and 
hence the reported rate. In future searches, it should be pos- 
sible to quantify the waveform uncertainties more precisely 
by performing simulations with amplitude corrected, spinning 
waveforms and by comparing the post-Newtonian waveforms 
to those obtained from numerical relativity simulations. 
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